Adhesion and electronic structure of graphene on hexagonal boron nitride substrates 
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We investigate the adsorption of graphene sheets on h-BN substrates by means of first-principles 
calculations in the framework of adiabatic connection fluctuation-dissipation theory in the random 
phase approximation. We obtain adhesion energies for different crystallographic stacking configura- 
tions and show that the interlayer bonding is due to long-range van der Waals forces. The interplay 
of elastic and adhesion energies is shown to lead to stacking disorder and moire structures. Band 
structure calculations reveal substrate induced mass terms in graphene which change their sign with 
the stacking configuration. The dispersion, absolute band gaps and the real space shape of the low 
energy electronic states in the moire structures are discussed. We find that the absolute band gaps 
in the moire structures are at least an order of magnitude smaller than the maximum local values 
of the mass term. Our results are in agreement with recent STM experiments. 



PACS numbers: 81.05.uc; 73.22.Pr; 71.10.-w; 71.20.-b 



I. INTRODUCTION 

The development of graphene-based nanoelectronic 
devices, such as high-speed transistors, calls for high 
electron mobilities. Currently, substrates beneath the 
graphene present an important source of disorder in- 
cluding corrugation effects of the graphenei^-4, charge 
trapsiiS"-^ and dangling bonds^*^'^^. A promising can- 
didate to become a new standard substrate material is 
hexagonal boron nitride (h-BN). This BN polymorph is 
remarkably similar to graphite: The alternating B and N 
atoms form two-dimensional layers of strong sp^ bonds 
within a honeycomb arrangement and a lattice constant 
which differs by less than 2% from that of graphene. The 
h-BN sheets are weakly bound by long-range adhesive 
forces at an equilibrium distance of S.SAiS,. The elec- 
tronic structure, however, exhibits clear differences: the 
chemically inequivalent sublattices make h-BN an insu- 
lator with a band gap of 6.0 eV^'^. Recently, the fabrica- 
tion of graphene devices on h-BN with highly improved 
electron mobilities and carrier inhomogeneities, reaching 
a quality comparable to suspended graphene has been 
reportedi^ii^. Thereby, the graphene was found to keep 
its zero band gap and to stack quasi randomly orientated 
on the h-BN substrate. 

In this paper, we analyze the adhesion behavior and 
the electronic structure of graphene on h-BN from first- 
principles. The paper is organized as follows: In section 
im we show that methods beyond standard density func- 
tional theory (DFT) are necessary to describe the weak 
non-local attraction between the h-BN and the graphene 
layers. We calculate adhesion energies using the ran- 
dom phase approximation (RPA) within the framework 
of the adiabatic connection fluctuation-dissipation theo- 
rem (ACFDT). On the basis of elasticity calculations, we 
discuss in section Hill mechanisms to release the stress re- 
sulting from the lattice mismatch and leading to the for- 
mation of moire superstructures. Section IIVI is devoted 



to the band structure and energy gaps of graphene on 
h-BN. From DFT band structure calculations we derive 
a low energy tight-binding description of graphene on h- 
BN and find mass terms which change their sign with the 
stacking configuration. This leads to an absolute gap in 
the moire structure which is at least an order of magni- 
tude smaller than the maximum local values of the mass 
term. The real space shape of the low energy states, 
particularly the issue of sublattice polarization and the 
occurrence of so-called snake states in regions where mass 
term changes its sign is discussed in section [V] Finally, 
conclusions and an outlook are given in section IVll 



II. ADHESION OF GRAPHENE ON h-BN 

DFT is a successful approach to describe ground state 
properties of solids. However, widely used semilocal ap- 
proximations for the exchange-correlation energy like the 
local density approximation (LDA) and generalized gra- 
dient approximation (GGA)^^'^^ do not take long-range 
correlations into account correctly. Thus, they fail to re- 
produce van der Waals attraction and the prediction of 
equilibrium geometries of van der Waals bound layered 
systems, such as graphite or h-BN, proves problematic 
with these methods^^ii^. A highly accurate means to 
describe van der Waals forces from first-principles is pro- 
vided by the random phase approximation to the corre- 
lation energjj22i^. Evaluated in the ACFDT framework, 
the RPA correlation energy reads2^ 
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where is the response function of the non-interacting 
Kohn-Sham (KS) system and v = YliKj 
electron-electron interaction. Together with the total 
KS Hartree-Fock energy, usually referred to as exact ex- 
change energy _Bexx^, the total ground-state energy of 
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the system reads as = E'exx + Ef^^. The ACFDT- 
RPA method has been proven accurate to describe bulk 
properties of solids such as lattice constants as well as 
adsorption energies^'^. For layered van der Waals bonded 
systems ACFDT-RPA yields a much more accurate de- 
scription of the structural properties than LDA and vdW 
DFT method&24'25. 

To simulate the graphene-h-BN system in this way, 
calculations were performed with the Vienna ab ini- 
tio simulation package (VASP)^^ using plane wave ba- 
sis sets within the projector-augmented wave (PAW) 
method^ii^. To this end, a unit cell (4 atoms) contain- 
ing graphene on an h-BN layer with 25A of vacuum 
above was constructed. Six stacking configurations were 
considered (Fig. [T^): starting from configuration I, the 
graphene sheet was translated downwards by half a B-N 
bond length in each step until the initial configuration 
was reached again. The lattice constant was chosen to 
2.49A, referring to the LDA optimized lattice constant 
of h-BN (a discussion of the lattice mismatch follows be- 
low), x^'^j Ef-^^ and E-^xx were evaluated with the 
LDA KS orbitals. In these computationally demanding 
simulations, convergence of the results was reached at a 
kinetic energy cut-off of 347 eV for the response function, 
a plane wave cut-off of 520 eV and a mesh of 7 x 7 x 1 
k-points. Additionally, standard LDA/GGA calculations 
were performed for the same geometries with a 24 x 24 x 1 
k-points grid and a kinetic energy cut-off of 500 eV. 

Fig. [5] shows the RPA total energies per unit cell 
(relative to the energy at large separation) for the dif- 
ferent stacking configurations as a function of the inter- 
layer spacing d. The curves show that, starting from the 
highest-energy configuration I, the lowest-energy config- 
uration V is approached step-wise. The configurations 
I, II and III exhibit total energy minima at interlayer 
spacings between 3.50A and 3.55A, whereat II and III 
are energetically virtually equivalent with minima of - 
65 meV; the minimum of configuration I is only slightly 
higher (— 62meV). Configuration V, where the carbon 
atoms sit on top of a boron atom and in the middle of 
the BN hexagon, is energetically most favorable with a 
minimum of -83 meV and an equilibrium layer distance 
of 3.35A. Energetically closest to V are IV (-71 meV) and 
VI (-70 meV), where the nitrogen atom is also not cov- 
ered by a C atom or a C-C bond. The curves exhibit no 
additional energy barrier for translation. For distances 
larger than 4A, all configurations become energetically 
indistinguishable. 

To get further insight to the nature of the attractive 
forces between graphene and h-BN, we analyze the decay 
of the RPA correlation energy Ef-^^ with the interlayer 
spacing d (Fig. [2] inset). We find E^^^ ~ d"*. This is 
clearly different from an exponential falloff as would be 
expected for local correlation effects as included in LDA 
or GGA. It rather indicates bonding of vdW type which 
yields power law decays^^ Ef'^^ ~ d^^ with p = 4 for 
2D insulating system&i^iiS. 

A comparison of the RPA calculations with the stan- 




Figure 1. (Color online) (a) Top view of the calculated stack- 
ing configurations for graphene on h-BN with the carbon 
atoms (yellow), boron (red) and nitrogen (light blue). Be- 
tween each neighboring images the graphene lattice is shifted 
downwards by half a B-N bond, (b) Moire structure with 
persisting lattice mismatch. A lattice mismatch of 1.8% cor- 
responds to a 55 x 55 moire unit cell. For clarity, a smaller 
moire unit cell (13x13) is shown, (c) Adhesion energy (-Eads) 
landscape in the moire pattern (color-coded), (d) The same 
for the local sublattice symmetry breaking A. 



dard LDA and GGA methods is given in Fig. [S] The 
LDA (see also '29') yields a qualitatively correct equilib- 
rium layer separation but underestimates the RPA ad- 
hesion energy by about 30%. The GGA, actually an 
improvement over the LDA in many casea^^^, exhibits 
an even more dramatic underbinding with a weak min- 
imum giving a negligible adhesion energy of less than 6 
meV/unit cell at an extremely high equilibrium distance. 
The LDA and GGA curves nearly coincide for distances 
larger than 4.5A and underestimate the van der Waals 
interaction also in this asymptotic region. The results 
clearly show that long-range correlations have to be taken 
into account for an accurate description of the layer at- 
traction. 



III. EXISTENCE OF MOIRE STRUCTURES 

We now turn to the discussion of the consequences of 
the lattice mismatch (1.8% in LDA and 1.9% in GGA) 
between the graphene and the h-BN and address the 
question whether stacking disorder or moire superstruc- 
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Figure 2. (Color online) Total RPA energies per unit cell for 
the stacking configurations I- VI as a function of the distance, 
d, between the graphene and h-BN layer. Inset: RPA correla- 
tion energy (relative to the correlation energy at large separa- 
tion) of configuration V as function of d~* between 2.9A and 
4.7A. 
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Figure 3. (Color online) Total energies per unit cell of config- 
uration V from RPA, LDA and CCA. 



tures should occur. To this end, we estimate the total 
energy difference of a structure with stacking according 
to the minimum energy configuration V in the entire sam- 
ple and a moire structure with persisting mismatch (Fig. 

In the case of persisting mismatch in the system, a 
moire pattern with a large unit cell (55x55 for 1.8% mis- 
match) is formed (Fig. [Ua). Here, out-of-plane corruga- 
tions resulting from interlayer spacings varying between 
3.35A (region V) and 3.55A (region I) are negligible, 
since their amplitude (0.2A) is small as compared to their 
wavelength (~135A). Hence, we focus on in-plane de- 
formations. It is visible from Fig. [T|d that our choice 
of stacking configurations I~VI simulated in RPA covers 
the sample uniformly and gives a sketch of the energy 
landscape (Fig. [!}:). Those parts of the moire pattern. 



where the nitrogen atoms are mainly beneath the cen- 
ter of the carbon rings (regions IV~VI) , are energetically 
more favorable than regions I-III. For the moire struc- 
ture we estimate the average adhesion energy per two 
carbon atoms by the average over the adhesion energies 
of the configurations I- VI. We obtain an average adhe- 
sion energy of 69 meV/(2 C-atoms). This is 14 meV/(2 
C-atoms) less than the adhesion energy of 83 meV/(2 
C-atoms) in configuration V. 

In the other case with the entire sample in configura- 
tion V, the lattice mismatch must be overcome and strain 
energy has to be brought up to force graphene and h-BN 
to have the same lattice constant. Now, two situations 
have to be dis ting uished. For graphene on a h-BN crys- 
tal as in Refs. [13! Isol. andlsi], the uppermost h-BN layer 
will likely keep its lattice constant at the bulk value. To 
stretch the graphene on the lattice constant of h-BN, 40 
meV/(2 C-atoms) of strain energy have to be overcome. 
Therefore, the strain energy overcompensates the adhe- 
sion energy gain of 14 meV/(2 C-atoms) by far when 
forcing the entire sample to configuration V. Hence, the 
lattice mismatch between the graphene and the h-BN will 
persist and strain will be released by realizing different 
stacking configurations as in the moire structure depicted 
in Fig. [T] This explains why multiple stacking coiifigura- 
tions have been realized in the experiment of Ref. [ij and 
also explains the recent observations of moire patterns in 
STM experiments^^. 

Second, one can conceive a situation where graphene 
is adsorbed on a free-standing monolayer of h-BN. In 
this case, the elastic properties of h-BN have to be ac- 
counted for since, then, a compression of h-BN is possi- 
ble. Our first-principles calculations (see App. yield 
two-dimensional Lame parameters and Young's moduli 
of A = 59N/m, fi = 125N/m, rh-BN,2D = 309N/m 
for single h-BN sheets. Compared to graphene, where 
1g,2D ~ 340N/m (yG,3D ~ 1.0TPa)32, the stiffness of h- 
BN is on the same order (about 10% smaller) and thus re- 
markably high. For graphene on free-standing h-BN, we 
find that a composition of stretched graphene and com- 
pressed h-BN is energetically most favorable and obtain 
a common optimized lattice constant of 2.467A (LDA) 
with the total strain energy being 18 meV/(2 C-atoms). 
Hence, this strain energy is very close to the adhesion 
energy gain of 14 meV/(2 C-atoms). This might lead to 
an interesting competition of these two energy contribu- 
tions and one might expect that systems with graphene 
on free-standing h-BN are highly sensitive to the experi- 
mental environment. 



IV. BAND STRUCTURE AND ENERGY GAPS 

We now investigate the band structure of graphene-h- 
BN hybrid structures and study the changes upon for- 
mation of a moire structure. To this end, we calculated 
the band structure of all geometries depicted in Fig. [T^ 
within the LDA. For configurations II, IV and VI we de- 
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I 


II 


III 


IV 


V 


VI 


A(meV) 
1 - t/t 


+57 



+7 
0.010 


-34 



-25 
-0.002 


-47 



-hl4 
-0.010 



Table I. Band gap A and ratio of the two inequivalent nearest- 
neighbor hopping parameters t/t with t = 2.45eV in struc- 
tures with broken trigonal symmetry. The average band gap 
is -4 meV. 



tect a small shift of the Dirac point in the hexagonal 
Brillouin zone away from the K to the M point (IV, VI) 
and in the opposite direction (II). Mapping the problem 
on a nearest-neighbor tight-binding model, we see that a 
description of graphene in configurations II, IV and VI 
requires two different hopping parameters, t and i, as 
the threefold symmetry of the graphene nearest neigh- 
bor bonds is broken — analogous to the case of uniax- 
ially strained graphene^. Fits of the TB model to the 
DFT results yield constant t = 2.45 eV in all regions and 
t 7^ i in regions II, IV and VI (Table H]) . In agreement 
with— we extract finite band gaps A in all regions 
varying between 7 meV and 57 meV. However, we find 
that these gaps have different signs (Table |l]), where we 
use the convention that a A > corresponds to states 
close to valence band maximum being entirely localized 
in sublattice B, while A < corresponds to states at 
valence band maximum being localized in sublattice A. 

In a moire structure like in Fig. [TJi this leads to a 
landscape of local sublattice symmetry breaking A with 
changing signs (Fig. [TJl, Table H]). We note that the local 
sublattice symmetry breaking does not necessarily lead 
to local spectral gaps in the LDOS. To gain understand- 
ing of the effect of the modulated gap landscape in the 
moire structure on graphene electrons, we consider the 
following tight binding model: 



a]ai 



bib, 



with t = 2.45 eV the nearest-neighbor hopping, a| (6|) 
the creation operators of an electron on sublattice A (B), 
and Ai the local mass term. The lattice vectors of the 
moire unit cell of size n x n are defined 
with di — a (1,0) and a2 — a (— 1/2, \/3/2) being the 
simple graphene unit cell vectors and a the graphene lat- 
tice constant. We denote local mass term by A^, where 
i = {l,in) describes the position within the the moire 
cell. Ai is periodic with the moire cell. Transforming the 
Ai to the reciprocal space, we find that the zeroth order 
Fourier component 



N ^ 



is given by the average of all local gaps in the moire 
cell. In addition to A^^^, the effect of the first order 
components of A^ with the smallest possible G ^ on 



the band structure is discussed in the following. To this 
end, we consider a sinusoidally modulated gap term 

Ai = Asin(27rZ/n + $i) -f B sin (27rm/n -I- $2) + C. 

(4) 

Here, B, C and $1,2 denote constants. Taking the lo- 
cal mass terms obtained from DFT in regions I- VI, realis- 
tic parameters A = 18.6 meV, B = 42.0 meV, $1 = 1.884 
and $2 = 1.531 can be obtained from a fit. 

First, we concentrate on the question whether the mod- 
ulated gap landscape opens an absolute band gap or not. 
In Fig. 2^, the two energy bands closest to the Fermi 
level for a 20x20 graphene supercell are depicted with a 
gap landscape as given by Eq. |4]and different amplitudes 
A, B . The average gap is set to zero (Ag^^ = (7 = 0). 
The green dashed line shows the case of a gap landscape 
with realistic amplitudes A, B as given above. We see 
that the bands nearly coincide with the bands obtained 
for entirely ungapped graphene (A^ = 0, red solid line) 
and most importantly, no band gap opens. This holds 
even for unrcalistically large values of the modulation 
amplitudes on the order of the hopping t (also Fig. |4^, 
blue and purple dashed lines). For large modulation am- 
plitudes, another important feature becomes visible: a 
renormalization of the Fermi velocity vp = j-^- For 
A, B being 100 times larger than the realistic values, up 
drops down by about 50% (purple dashed line). For a 
realistic gap landscape, though, this effect is too small 
to be detected in experiments — in contrast, e.g., to the 
case of twisted bilayer graphene^. 

So, near the Fermi level, the energy bands of graphene 
are only weakly affected by a modulated gap landscape 
with realistic parameters and zero average gap - no band 
gap opens and the amplitudes are too low to renormal- 
ize the Fermi velocity measurably. At the Brillouin zone 
boundary (Fig. HJd), however, a difference to the per- 
fect isolated graphene becomes apparent: here, mini- 
gaps open similar to the case of graphene on Ir(lll) 
moires'^^''^®. 

Now we discuss the scenario of a non-zero average gap, 
i.e., the case where the zeroth Fourier component is non- 
vanishing (Ag^p — C ^ 0). Fig. [S] shows the bands 
near the Fermi level of a realistic gap landscape with 
C — 5 meV (green dashed line) and C = 15 mcV (blue 
dashed line). Here, an absolute band gap on the order 
of C opens that remains stable upon adding Ag_^p terms 
of realistic amplitudes. Similar as for the Fermi velocity 
discussed in Fig. |4] (a) , our calculations showed that the 
band gap reduces measurably (but not entirely closes) 
when the modulation amplitude is increased by orders of 
magnitude. However, in the realistic scenario of graphene 
on h-BN, Fig. [5] clearly shows that the only quantity de- 
termining the absolute band gap is the zeroth Fourier 
(■3^ component A^^^p. Hence, the average gap Ag^^ corre- 
sponds to the absolute spectral gap in the moire struc- 
ture, while amplitude and periodicity of spatially oscil- 
lating contributions to A renormalize the Fermi velocity. 
For the structure of Fig. [Hd we find A w —4 meV. There- 
fore, we expect an absolute gap which is at least an order 
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Figure 4. (color online) (a) Bands close to the Fermi level of a 20x20 moire cell of entirely ungapped graphene (red line) and 
graphene with sinusoidally modulated gap (average gap ^q^q = 0) with realistic amplitudes (green dashed line) and increased 
amplitudes A, B (blue/purple dashed lines). No band gap opens. For visualization purposes, the Dirac point has been shifted 
back to the K point, (b) Bands below the Fermi level at the Brillouin zone boundary (folded to the M point) for ungapped 
graphene (red line), realistic amplitudes (green dashed line) and increased amplitudes (blue dashed line): minigaps open for 
finite amplitudes. 




0.346 



Figure 5. (color online) Bands close to the Fermi level of 
a 20x20 moire cell with a sinusoidally modulated gap with 
realistic amplitudes and non-zero average gap Ag^Q = C. 
Band gaps as large as C open. 



of magnitude smaller than the maximum local values of 
|A|. This is well in line with the absence of a gap being 
reported in transport^ and STM experiments^. Our TB 
simulations further show that velocity renormalizations 
are below 2% for the moire structure of Fig. [T|d. 



V. REAL SPACE SHAPE OF LOW ENERGY 
STATES 

To understand how spatially modulated gap terms 
(Ag_^p) affect the graphene electrons and how they mani- 
fest, for instance, in local probe experiments, we visualize 
the states close to the Fermi level in real space. The Figs. 



[6l [7] and [8] illustrate a 20x20 graphene supercell with sub- 
lattice A atoms as dots in black and sublattice B atoms 
as dots in red color. Here, the size of the dots illustrates 
the contribution of each atom to the states below the 
Fermi level (Figs. EHH^) and above (Figs. HHHb) in close 
proximity to the Dirac point. 

Fig. [5] shows the case of realistic modulation ampli- 
tudes and vanishing average gap A^^^ = 0. In this 
case, the amplitudes of low energy states are equally 
distributed over both sublattices throughout the entire 
moire cell as in ungapped graphene. There is no clear 
enhancement or decrease of probability density in any 
region of the moire cell. However, increasing the gap 
modulation amplitude by a factor of 25 (Fig. [7]) induces 
a localization of the low energy states in regions where 
the local gap |Ai| is small — so-called snake states occur, 
but no absolute band gap opens. Apparently, the states 
are equally localized in both sublattices. 

Whether or not snake states occur depends on the ratio 
of the modulation amplitudes A, B to the energy En ~ 
27ThvF /(na) related to the moire periodicity na. 20x20 
and 50x50 moire cells lead to £"„ « 0.7 eV and 0.3 eV, 
respectively. In the case of large modulation amplitudes 
A,B > En (corresponing to Fig. [T]), the wave functions 
near the Dirac point clearly have the shape of snake states 
— in contrast to the case of A,B <C En, where the low 
energy LDOS is almost homogeneous in the entire moire 
cell (Fig. [6]). The situation of Fig. [6] is well in line 
with no gap or LDOS inhomogeneities being detected in 
the local probe experiments of Ref. [s^- It is, however, 
important to note that a moire periodicity of 100a which 
might be reached by external strain or twisting already 
corresponds to an intermediate case of A, B, and £"„ 
being on the same order magnitude. 

The shape of the low energy states again changes for 
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Figure 6. (color online) 20x20 graphene supercell (red: sublattice A, black: sublattice B) with sinusoidally modulated gap 
terms with realistic amplitudes of A = 18.6 meV, B = 42.0 meV and vanishing average gap (Ag^p = 0). The size of the dots 
depicts the contribution of each atom to states close to the Dirac point in an infinitesimal energy window around the valence 
band maximum (a) and the conduction band minimum (b). 
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(b) 



Figure 7. (color online) The same as Fig. [6]for a 20x20 moire cell with a sinusoidally modulated gap with amplitudes A — 25A, 
B — 25B, and zero average gap (A^^^ — 0). So-called snake states occur. 



a system with a non-zero average gap A^^^ fa A, B 
on the order of the modulation amplitudes. In Fig. [5] 
the case of a realistically modulated gap landscape and 
Ag^P — +50 meV is depicted: the states are homoge- 
neously distributed in space fully but sublattice polar- 
ized - with the state below Fermi level almost entirely 
localized in sublattice B (Fig. [5] (a)) and vice versa 
(Fig. [8] (b)). Our DFT calculations yield the case of 
a much smaller average gap on the order of few meV 
(A^^Q <^ A,B). In that case, two situations have to be 
distinguished: First, in the case of A^ „ A, B En 
sublattice polarized states as in Fig. [8] should be de- 
tectable in STM experiments, if the LDOS is measured 
inside an energies range E < A^^p which is sufficiently 
close around the Dirac point. Otherwise low energy 
states without any sublattice polarization as shown in 
Fig. |6]will be detected. This situation would be again in 
line with the STM experiments of Ref. [s^. 

Differently, in the case of ^q^q ^ En ^ A, B, there 
is generally no the sublattice polarization detectable in 
the low energy LDOS but snake states similar to Fig. [7] 
occur. 



VI. CONCLUSIONS 

In summary, we have calculated accurate adhesion en- 
ergies for graphene-h-BN systems by means of ACFDT- 
RPA. A comparison of the strain energies with the adhe- 
sion energy differences suggests that a lattice mismatch 
between h-BN and graphene persists in experiments like 
Refs. 13 [sO, and [Sll This explains the experimental 
observation of different stackings, moire patterns, and 
stacking disorder. Our band structure calculations show 
that the gap landscape in the moire structure exhibits 
mass terms with changing sign and a small average gap. 

The interplay of constant and spatially oscillating gap 
terms is decisive for determining whether or not phe- 
nomena like Anderson localization can occur— — partic- 
ularly in experiments where graphene is very close to the 
charge neutrality point^. Gaps with spatially changing 
sign also control the transport properties of systems like 
(Hg,Cd)Te quantum well structures^ which can be tuned 
into a topological insulator. While we find that struc- 
tures like those in Refs. [ij, [s^, and [sll are likely not in 
a "topological insulator regime" , where charge transport 
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Figure 8. (color online) The same as Fig. [6] for a 20x20 cell of sinusoidally gapped graphene with realistic amplitudes A, B 
and a finite average gap (here: Ag^^ = 50 meV). 



would occur through protected edge states, it remains 
to be seen whether this might be reaUzed in structures 
like externally strained graphene on free-standing h-BN, 
where considerably larger moire periodicities might be 
realized. 
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Appendix A: Elastic properties of h-BN sheets 

Here, we now discuss the calculations of the elastic 
constants of h-BN (the elastic properties of graphene 
have been widely investigated in experimen t'^^'^^ and 
theory^"). The isotropic Young's modulus and Lame 
parameters of single h-BN sheets were obtained from 



DFT calculations. These elastic constants are deter- 
mined by strong in-plane chemical bonds and well de- 
scribed within the LDA/GGAii. For single h-BN sheets, 
the two-dimensional Young's modulus is defined by 



l2D = 



1 d-'E, 



(Al) 



where e is the axial strain, Es the total strain energy 
and ^0 the equilibrium surface. The strain energies of 
a h-BN sheet were evaluated with uniaxial strains be- 
tween -8% (compression) and 8% (tension). We obtain 
the two-dimensional Lame parameters and Young's mod- 
uh of A = 59N/m, ^ = 125N/m and yh-BN,2D = 309N/m 
within LDA and A = 54N/m, = 123N/m, Yh-BN,2D = 
300N/m within GGA. The ih-BN,2D correspond to three- 
dimensional Young's moduh of yh-BN,3D = 0.94TPa 
(LDA) and lh_BN,3D = 0.91TPa (GGA), assuming an 
interlayer separation of B.sAi^. Our full potential PAW 
calculations yield about 10% higher elastic constants 
i^-BN,2D and ^ than calculations using Gaussian basis 
sets^. The results are in agreement with the experi- 
ment, where yh-BN.2D,oxp ~ 220 — 510N/m was obtained 
for few-layer h-BN'''^. 
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